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Among the various numerical techniques to study the physics of strongly correlated quantum 
many-body systems, the self-energy functional approach (SFA) has become increasingly important. 
In its previous form, however, SFA is not applicable to Bose-Einstein condensation or superfluidity. 
In this paper we show how to overcome this shortcoming. To this end we identify an appropriate 
quantity, which we term D, that represents the correlation correction of the condensate order pa- 
rameter, as it does the self-energy for the Green's function. An appropriate functional is derived, 
which is stationary at the exact physical realizations of D and of the self-energy. Its derivation 
is based on a functional-integral representation of the grand potential followed by an appropriate 
sequence of Legendre transformations. The approach is not perturbative and therefore applicable 
to a wide range of models with local interactions. We show that the variational cluster approach 
based on the extended self-energy functional is equivalent to the "pseudoparticle" approach in- 
troduced in Phys. Rev. B, 83, 134507 (2011) We present results for the superfluid density in the 



two-dimensional Bose-Hubbard model, which show a remarkable agreement with those of Quantum- 
Monte-Carlo calculations. 

PACS numbers: 64.70.Tg, 67.85.Dc, 03.75.Kk 



I. INTRODUCTION 



Seminal experiments with ultracold gases of atoms 
trapped in optical lattices shed new light on strongly- 
correlated many body systemsi^"— In these experiments 
specific lattice Hamiltonians can be engineered and inves- 
tigated with a remarkable high level of control, making 
quantum mechanical interference effects observable on a 
macroscopic scale. Most important as well as fundamen- 
tal is the quantum phase transition of strongly corre- 
lated lattice bosons from the localized Mott phase to the 
delocalized superfluid phase. In the superfluid phase a 
macroscopic fraction of the particles condenses into one 
quantum mechanical state, thus, forming a Bose-Einstein 
condensate, where the number of particles in the con- 
densate is not necessarily equal to the number of super- 
fluid particles. In experiments with ultracold gases of 
atoms trapped in optical lattices, the condensate den- 
sity can be extracted from time-of-flight imageSf^ which 
are related to the momentum distribution of the conflned 
particles. Importantly, the flnite expansion time of the 
particle cloud has to be taken into account when draw- 
ing the connection between these time-of-flight images 
and the true momentum distribution4i"— However, it is 
probably even more challenging to measure the super- 
fluid density itself, as it is not a ground state property 
but rather related to the response of the system to a 
phase twisting field Interestingly, only very recently for 
Bose gases without the periodic lattice potential an opti- 
cal method has been proposed to extract the superfluid 
density. This experiment creates a vector potential, 
that imposes angular momentum on normal fluid parti- 
cles, while superfluid particles stay at rest<^ 

In a previous work, we extended the variational cluster 



approach (VGA) , which is capable to deal with strongly- 
correlated many body systems without broken symme- 
try to the superfluid phase of lattice bosons Originally, 
VGA has been formulated for the normal Mott phase of 
lattice bosons in Ref. [Ill within the so-called self-energy 
functional approach (SFA), which was previously intro- 
duced for interacting fermionic systems.— Our exten- 
sion to the superfluid phase in Ref. [13 follows a different 
path, and is based on the so-called "pseudoparticle" for- 
malism. Within this approach we obtained the expres- 
sions for the superfluid order parameter, the anomalous 
Green's function, and the grand potential, which is the 
starting point for the variational principle, see Eq. (1), 
(33), and (2) in this reference. 

It should be pointed out that, while the pseudoparticle 
formalism is equivalent to VGA in the normal phase of 
both bosonioiS and fermionioi^ systems, it lacks the rig- 
orous theroretical framework provided by SFA. In partic- 
ular, there is no genuine variational principle explaining 
why one should look for a saddle point in the grand po- 
tential. The goal of the present paper is to put the results 
obtained within the pseudoparticle approach into a rig- 
orous framework by developing an extended self-energy 
functional approach, which is capable to deal with the 
bosonic superfluid phase. 

From the present work it will become clear, that this 
extension is not straightforward, as it involves a precise 
sequence of Legendre transformations with suitably cho- 
sen variables. In the search for the appropriate set of 
transformations the knowledge of the final results pro- 
vided by pseudoparticle formalism proves to be useful. 
This fact emphasizes the advantage of the heuristic, yet 
straightforward, pseudoparticle approach to formulate 
extensions of VGA 

The extended SFA formulated in the present paper 
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yields the same expressions for the superfluid order pa- 
rameter, for the Nambu Green's function, and for the 
grand potential, as obtained from the pseudoparticle ap- 
proach. While this might not seem to be surprising, since 
we were guided by the very results of the pseudoparticle 
approach, we argue below, that our SFA extension pre- 
sented here is unambiguous. The most important step 
in this SFA extension is to find a quantity, which we 
call U, which is the companion of the self-energy in the 
superfluid phase. Correspondingly, one has to find an 
appropriate universal functional of this quantity and of 
the self-energy, which generates the superfluid order pa- 
rameter and the Green's function. 

As an application, we present an evaluation of the 
superfluid density within this extended VGA, by the 
usual method of introducing a phase twisting field, 
which is equivalent to the helicity modulua^^ and to 
winding numbers in quantum Monte Carlo (QMC) 
algorithmsii^iii We evaluate the superfluid density for 
the two-dimensional Bose-Hubbard (BH) modeli^ and 
find a remarkably excellent agreement with QMC results. 

This article is organized as follows. In Sec.|ll]we extend 
SFA to the superfluid phase and obtain the corresponding 
extended self-energy functional, along with the appropri- 
ate variable describing superfluidity. The evaluation of 
the superfluid density within this extended VGA is pre- 
sented in Sec. IIIII and applied to the BH model in two 
dimensions. The VGA results arc compared with unbi- 
ased QMC results showing excellent agreement. Finally, 
we conclude and summarize our findings in Sec. IIVI 



II. SELF-ENERGY FUNCTIONAL APPROACH 

Let us recall the key idea of SFA due to M. PotthoffJ^ 
The starting point is an appropriate functional 

f^p, G^\Hu] = J-p, Hu] + Go 1] , (1) 

which consists of a functional of the self-energy, the 
Legendre transform of the Luttinger-Ward functional, 
which is universal in the sense that it depends on the 
interaction part {Hjj) of the Hamiltonian but not on the 
single particle part. The latter enters via the free Green's 
function Gq^ in the second functional, which is explicitly 
known 

Go = trln(S]-Go"i) 

The functional r2[S, Gq"^, has three key features, 
which arc crucial for VGA. 

a) The non-universal part £ enters additively in form 
of a known functional and the many-body aspects 
are described by a universal functional independent 
of the single particle Hamiltonian, or, cquivalently, 
independent of Gq^^. 



b) The self-energy of the physical system, character- 
ized by Hu and Gp ^ is a stationary point of the 
functional 17 with respect to S. 

c) The value of at the stationary point is equal to 
the thermodynamic grand potential. 

Given these properties, one can construct a paramet- 
ric family of Hamilton operators based on the same in- 
teraction part (reference systems) , for which the thermo- 
dynamic grand potential, the Green's function and the 
self-energy can be determined exactly. This allows to de- 
termine the exact self-energy functional for self-energies 
accessible by the reference systems. In this very sub- 
space, the self-energy functional in Eq. ([T]) for the phys- 
ical system is replaced by that of the reference system. 
The stationarity condition in turn allows to determine 
the Green's function and self-energy of the physical sys- 
tem. 

Our goal is to generalize this approach to the super- 
fluid phase as well. Besides the self-energy, which is the 
interaction correction of the inverse Green's function, we 
need the corresponding companion that describes the in- 
teraction correction to the order parameter, which we call 
D. 

Once the appropriate form of D has been determined, 
we need a functional 

(ls[^,D,F,G^\Hu]^P[^,D,Hu]+£[^,D,G^\F] , 
in the self-energy E and D with the following features. 

a) J- is again a universal functional, now in E and D. 
The non-universal part E is explicitly known and 
carries the dependence on Gq ^ and the symmetry 
breaking source-field F . 

b) The functional is again stationary at the exact self- 
energy S and the exact D of the physical system, 
characterized by Hu^ Gq^ and F. 

c) The value of Qg at the stationary point is equal to 
the thermodynamic grand potential. 

The sought- for functional fig, to be derived in this sec- 
tion, will turn out to be (see below for a definition of the 
quantities) 

2/?f7,[i], d,Go\f]^t[j:,d] + £[j:, d,g^\f] (2) 

D,G^\F]=pTv ln[(Go 1 - E) Goo] 

+ {D^F){G^^-^)-\D^F). 

(3) 

In the normal phase, it is identical to the functional 
introduced by Potthoff. The additional factor 2 is due to 
the Nambu Green's functions. Moreover, the expression 
for the grand potential obtained with the help of a so- 
called reference system, see Eq. [3D] below, is identical to 
the one obtained within the pseudoparticle approachiiS 
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A. Derivation of the grand potential functional 

We start out from the partition function Z of a bosonic 
many-body system, which in a functional integral repre- 
sentation reads 



Z = VA 



(4) 



where S is the action, which in general can be written 



S^-^ I dr I dr'A{r')G-\T',r)A{T) 



dr [Fir) Air) - HuiA{r))] 



(5) 



In view of treating the supcrfluid phase we have adopted 
a Nambu notation in which the bosonic fields are ex- 
pressed in a vector representation 



Air) 



a-Nir) 
ai(r) 



V aNir) J 



(6) 



The indices 1 through N denote the corresponding single- 
particle orbitals (for example, lattice sites) where the bo- 
son operators act, and aj;(r) idi{r)) are the fields asso- 
ciated with the annihilation (creation) of a boson in the 
orbital i. The adjoint field is defined as 

Air) EE (ai(T), ■ • • ,aNir),aiir), ■ ■ ■ ,aAr(r)) . (7) 

It can be expressed in terms of Air) with the help of the 
matrix T, which exchanges the first entries of a vector 
with the last TV ones: 



Air) = Air)^T . 



(8) 



The operator T has the properties = 1 , and T = . 
The action in Eq. (O also contains the source fields 

F= (/i,.../^,/i,... ,/Ar) and F = TF^, 

which are zero for the physical system of interest, the bo- 
son interaction described by Hu, as well as the 27V x 2N 
noninteracting Green's function matrix Go(t', r). Eq. (|4]) 
with Eq. JS]) defines the corresponding grand potential as 
a functional of Gn^ and F 



1 



ns[Go\F] = --\nZ 



(9) 



where /? is the inverse temperature. Here and in the fol- 
lowing, we mark functionals with a hat " " " ^ and omit 



their arguments whenever they are obvious. The nonin- 
teracting Green's function has the matrix structure (see 
App.Ell) 



G^\r',r) = -5ir-r')(^ 



dr 











-dr 



(10) 



where t is the single-particle Hamiltonian matrix. 

In the following, we carry out a sequence of Legcndre 

transformations starting from fig , ultimately leading to a 
universal functional J^[S, D] of the self-energy E and of a 
suitable quantity D defined in Eq. (|21ap . The functional 
T is the generalization of the self-energy functionalii"— 
to the superfluid phase, where a nonvanishing expecta- 
tion value Air)={Air)) of the boson operators A exists. 
The functional has the properties, see Eq. that its 
functional derivatives with respect to S and D yield the 
disconnected Green's function, and the expectation value 
A, respectively. This procedure is inspired by Ref. [13 and 
extends that approach to the treatment of the superfluid 
phase. 

We first determine the conjugate variables to Gq^ and 
to the source fields F. The functional derivative of 
fls with respect to the noninteracting Green's function 
yieldai^ (sec App. lAip 



2/3 



VAx 



^Goj,(T',r) Z<5Goj,(r',r) 

xe^v{\jdr j dr'Aiir)G-^lir,r')A,ir') 

dr[Fiif) A^ ir)~Huia, a)]} 

1 



= -jj PA^,(r')A,(r)exp[-5] 

= Gdisc,ij(''", r') . 

Here Gdisc,^ (t, t') is the disconnected interacting time- 
ordered Green's function. Along with the definition of 
the connected Green's function G[Gq^,F] we obtain 



2/3 



6ns[G^\F] 
SGn' 



Gd 



G-AA 



(11a) 



For the functional derivative with respect to F we obtain 
similarly 



2/3^^fi^ = -2i[G-,F]. 



(lib) 



The two functionals G[Go \ F] and A[G^'^,F] defined in 
Eq. (fTTj) provide the exact Green's function G and order 
parameter A for a given noninteracting Green's function 
Gq^ and source field F of the system. The first step 
toward the universal functional consists in a Legendre 
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transformation replacing the variables F with A. To this 
end, we invert^ the relation Eq. pT|) making F a func- 
tional F[Gq^,A] and introduce 



(12) 



where, as usually in Legcndre transformations, the func- 
tional dependence on F has been eliminated in favor of 
A by using the inverse function. It is straightforward to 
show that the corresponding functional derivatives give 



52A 



F[G^\A] 



6Gn 



Gdisc[Go ^,A] 



Next, we modify the functional in the following way 

i[Go\A] = E + AG^'A, (13) 

such that we obtain the connected Green's function from 
the functional derivative with respect to Gq ^- In total 
we have 



SE 
62A 



= f[Go\a] + AGo 



SE 



SGo 



— = Gdisc +AA = G[Go\A] 



(14) 
(15) 



The second step is a Legendre transformation replacing 
the variable Gn^ with G 



n[G, A] = S-/3Tr(GGo 



1) 



(16) 



where we have expressed Gq ^ as a functional of G and 
A, by inverting Eq. ()15p .— i2i We subtract an "infinite" 
constant /? Tr 1 in order to keep II[G, A] finite. The func- 
tional derivatives of the new functional are 



— -F + AG, , — 



Gn 



Now, we modify the functional such that we get the self- 
energy from the functional derivative (see App. lA 1 cp 



n[G, A] =n + ^Trln(G/Goo) 



This gives 



S2A 



Yj — G Gn 



(17) 

(18a) 
(18b) 



Here we have used the Dyson equation as defining equa- 
tion for the self-energy. Furthermore, we carry out a 
third Legendre transformation replacing G with E in the 
usual way. Thus we introduce 



P[E, A] = n + /3TrEG 



(19) 



with the properties 
SP 



52A 



F 



-AG^' 



SP 



G 



We modify this functional once more so that its derivative 
yields a new function D, which will be the companion of 
the self-energy in our extended self-energy approach 



A] = P- AT.A . 
The functional derivatives yield 



(20) 



AP - _ _ ^ _ ± 

— =F + AG^^-AJ: = F + AG-'^=D, (21a) 
S2A ' y ' 

SP 



<5S 



^G-AA^G 



disc 



(21b) 



Before proceeding, let us discuss the meaning of the 
function D introduced in Eq. (|21ap . When extending 
SFA to the superfiuid phase one is looking for a quan- 
tity, which is related to the condensed order parameter 
and which plays a similar role as the self-energy, in that it 
describes the deviation between the interacting and non- 
interacting case. Thus, this quantity should vanish in the 
noninteracting case {Hu = 0). The reason is that SFA 
will eventually amount to an approximation for E and D, 
and we require this approximation to become exact for 
Hij — 0. Finally, D must obviously vanish in the normal 
phase. The expression in Eq. (|21ap has precisely these 
features, since Aq = —FGq, which is straightforwardly 
determined from the Gaussian integral for Hu = in 
Eq. ([5|) . Interestingly, the pseudoparticle approach, pre- 
sented in Ref. [lol . and which is based on an intuitive, yet 
heuristic approximation, provides the same form of D as 
given in Eq. pia^ . 

The final Legendre transformation replacing A with D 
yields the desired functional of the self-energy and D. It 
represents the generalization of the self-energy functional 
(F[E] of Refs. [ill and [Til) to the superfiuid phase 



J^[E, D]=P - 2DA 
and has the properties 



(22) 



S_T_ 
W 



-2A[j:,d] , 



ST 



Gdisc[E,Z?]. (23) 



Similarly to F[Y,] from Refs. [T3 and [ill J- is (for fixed 
Hu) a universal functional of E and D only, from which 
the disconnected Green's function and the order param- 
eter are obtained by functional derivative, see Eq. ((23)) . 

Given E and D we can compute by Eq. ([^ the corre- 
sponding values for A and Gdisc- On the other hand, for 
a specific physical system, uniquely defined by Gq^ , F 
and Hu, the definitions of the self-energy E, Eq. (jl8bp . 
and the modified order parameter D, Eq. (|21a[) . provide 
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another set of equations, which uniquely fix E and D via 
the equations 

Gdisc[S,^] = (Goi-S)-i+ (24a) 
(G„-i - ^)-\D F){D - F)(Go 1 - S)-i , 



and 



^2A[Y., D] = -2{D - F){Go^ - S)" 



(24b) 



As for the (original) self-energy functional approach, we 
seek now a functional, which becomes stationary at the 
exact E and D for specific and and which consists 
of the universal functional J- plus a non-universal explicit 
functional of the form 

2l3n, [E, D,Go\ F, Hu] - -F[E, D, Hu] + 5[E, D, F] 

In order to yield the correct stationary point, the func- 
tional £ has to fulfill according to Eq. ([^ the equations 



(25a) 



(Go 1 - ^)-\D F){D - i?)(Go-i - E)-i , 



|^ = 2(^-^)(Go-i-E)-i 



(25b) 



With these ingredients we can now express the sought- 
for functional Cls as 

2/?f7,[E, Z?,Go\^^] = J-[E,i?]+/3Trln[(Goi -E) Goo] 

+ {D-F){G^'-^)-\D-F), 

(26) 



which obviously fulfills Eq. ([25|) . It remains to show that, 
whenever evaluated at the exact E and D the functional 
Vts corresponds, possibly apart from a constant, to the 
thermodynamic grand potential fig of the system. To 
this end we add up all the terms used to construct the 
functional. At the exact values of E and D we have 



-1 





2/3f^-L.act = + 2^-^ + -^Go M - /3 Tr (GG, 

+/3 Tr In (G/Goo ) + /? Tr EG - ^E^ - 2b A 
-/3Trln(G/Goo) +^G-U 
= 2/3f2g - 2{D - F)A + 2AG-^A 

= 2phs 

We can now proceed as in Refs. [l^ and[22l and construct a 
reference system, which can be solved (almost) exactlyi^ 
The reference system is described by a Hamiltonian H', 
which shares the same interaction Hu as the physical 
system, but consists of different noninteracting Green's 
function Gq and source fields F' . The point is the follow- 
ing: Due to the fact that is a universal functional, it 



cancels out from the difference between fJ^ for the phys- 
ical and the reference system, with the same values of E 
and D. In particular, this gives 

2/317,[E, D, Go \ F] - 2/3i],[E, D, G^"\ F'] 
= /3TVln ((Go-^ - E)Goo) - /3Trln ((G;,-! - E)Goo) 
+ {D~F){G,'-J:r\D-F) 

-{D-F'){G'^-'-^)-\D-F'), (27) 

which allows to evaluate the functional VLs exactly for the 
physical system as well, however, in a restricted subspace 
of E and D, representable by the parametric family of 
reference systems. By construction, the optimal values 
for E and D for the physical system are those of the 
reference system for the set of optimal variational pa- 
rameters. 

The variational procedure then follows and generalizes 
Rcf. [13 First a class of exactly solvable reference systems 
H' with the same interaction as the physical system char- 
acterized by a continuum of single-particle parameters 
t' and source fields F' is identified. In VGA this class 
is obtained by dividing the original lattice into discon- 
nected clusters with varying single-particle energies and 
hopping strengths. A larger subspace can be reached by 
adding bath sitesJ^ Then the (connected) Green's func- 
tion G', the order parameter A! , and the grand potential 
il'g of the reference system is evaluated. With the help of 
Dyson's equation Eq. (jl8bp the self-energy E', and with 
the help of Eq. (|21ap D' is determined. By varying t' and 
F' the subspace of self-energies and Ds is spanned, which 
is accessible to the reference system and to which these 
objects for the physical system are restricted. Within 
this subspace the functional 1)^ can be evaluated exactly 
for arbitrary Go and F of the physical system. For the 
relevant case F = we obtainSi from Eq. (P7)) 



2l3ns = 2/317; 4-/3 Tr In (-(Go"^ - E')) 

-/3Trln(-(G;,-i-E'))+^(G 
- A!G'-^A , 



-1 





(28) 



which is now a function of t' and F'. The infinite physi- 
cal system can break the symmetry spontaneously, while 
in the reference systems of disconnected finite clusters, a 
non- vanishing order parameter can only be achieved by 
an additional source field F'. This explains, why a finite 
F' is required although F = in the physical system. 
The SFA approximation consists in finding a stationary 
point of within this subspace of self-energies and D-s. 
This corresponds, quite generally, to finding a stationary 
point with respect to t' and F' of Eq. (|28p . i.e. to the 
equations 



OF' 







(29) 



Here, we have replaced ilg with fl = fig — ^ tr t which dif- 
fers just by a t'- and F'-independent constant and thus 
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does not change the saddle-point equations. The quan- 
tity ri is the grand potential obtained from the normal- 
ordered Hamiltonian (see App. lA 3[) . We also introduce 
the grand-potential of the normal-ordered reference sys- 
tem Jl' = — i trt'. This term is also present in the 
pseudoparticle approach^i^ where its origin is easily seen. 
Moreover, for r-independent fields and Hamiltonian, the 
expectation values A{t) are r-indepcndent as well, and 
the Green's functions depend on the time difference only. 
In this way, we can rewrite Eq. 



n = n' - itr(t - t') - iTrln(-G) + ^Ti-\n{-G') 
+ ^AG-' (w„ = 0)^ - i^'G'-i (w„ ^ 0)A' , (30) 

where G(a;„) = J d t G(t, 0)6'"^"" is the Green's func- 
tion in Matsubara space. The expression for f2 given 
in Eq. (|30p is our main result. As can be seen, this ex- 
pression is the same as Eq. (1) in Ref. [l^, except for 
a different normalization factor, which is the number of 
clusters Nc- Notice that Nch in Ref. [13 is equal to t — t' 
in the present paper. We thus proved that the result ob- 
tained within the pseudoparticle approach in Ref. [l^ can 
be equivalently obtained within a more rigorous "gener- 
alized" self-energy functional approach. While the pseu- 
doparticle approach is quite intuitive, the present self- 
energy approach provides a rigorous variational principle, 
explaining why the grand-potential 57 has to be optimized 
with respect to the cluster parameters t' and F' . In ad- 
dition, as in SEA for the normal phase, it suggests more 
general approximations in which bath sites are used to 
enlarge the space of possible self-energies J^i^i^^ 



III. SUPERFLUID DENSITY 

In this section we discuss the evaluation of the super- 
fluid density ps within our extended SEA/ VGA theory 
and present results for the two-dimensional BH model. 

The superfluid density is related to the response of the 
system to a phase-twisting fieldj ^^i^^ leading to twisted 
boundary conditions (BG) in one spatial direction, which 
we choose to be the e^j-direction, and periodic BG in the 
others. The many-body wave function \^) has to obey 
these BG and thus 

f(iV,e,)|vI') =e'«|vl/) , 

where the operator T(r) translates the particles by the 
vector r, is the lattice extension in e^-dircction, and 
is the phase twist applied to to the system. The twisted 
BG can be mapped by a unitary transformation onto 
the lattice Hamiltonian, leading to complex-valued hop- 
ping integrals j^i^iS^ The resulting Hamiltonian can be 
interpreted as a cylinder rolled up along the x-direction, 
which is threaded by an effective magnetic field with total 
fiux 8. When a particle is translated by Nx in the e^,- 
direction a phase exp[— i6] is picked upi^ Due to gauge 



invariance, one is free to choose where the phase is col- 
lected when the particle propagates across the lattice. 
The usual choice is that each hopping process in the e^, 
direction, i. e., from site r' = (r^, — 1, ry) to r = (r^,, Vy), 
is multiplied by a phase factor exp[— lA], where the asso- 
ciated vector potential is 



A = Q/Nx 



(31) 



When choosing the phase in that way, the reference 
system H' also depends on the vector potential A and 
the intra-cluster hopping terms become complex-valued 
along the e3;-direction. Eor a Hamiltonian with nearest- 
neighbor hopping i, the superfluid density is determined 
from^ 



1 1 17e 

iNxNy dA^ 



(32) 



where Ny is the total number of lattice sites of the 
physical system, and i^e is the grand potential of the 
physical system, subject to a phase twist 0, as discussed 
above. Plugging in the vector potential of Eq. (PT|) yields 



tNy 992 



(33) 



In practice, the grand potential fie is evaluated at the 
stationary point of Eq. pop , and is determined self- 
consistently for several values of &. Erom this data the 
curvature of fie with respect to O is extracted from a 
fit. Using the curvature, the superfluid density is evalu- 
ated according to Eq. ([33|) . Note that a finite cluster is 
embedded in an essentially infinitely large system and 
thus the limits are taken in the correct order to obtain 
the superfluid density.— 

In the following, we apply this procedure to the two- 
dimensional BH modeli^ 



H 



(i.j) 



^ n, [hi 



where a\ (a,) creates (destroys) a bosonic particle on site 
I, and fii = a\ is the occupation number operator. The 
hopping integrals tij are nonzero for nearest neighbors 
only, as indicated by the the angle brackets. Specifically, 
tij ~ —t for hopping processes along the e^-direction 



and ti 



-tcyi]i[i A{vi — fjOe^r] for hopping processes 



along the e^j-direction. The chemical potential, termed 
/i, controls the particle number and U is the repulsive 
on-site interaction, which subsequently will be used as 
unit of energy. The reference system H' consists of a 
cluster decomposition of the physical system H plus a 
U{\) symmetry breaking source term 

R (q,0) a 

- /i' ^ n„,R - ^ (a|^ /* a^ ^) 
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FIG. 1: (Color online) Superfluid density pa (a) evaluated for constant hopping strength t/U — 0.02 as a function of the 
chemical potential p/U. VGA results for reference systems of size L = lxl,L = 2xl, and L = 2 x 2 and for essentially 
infinitely large physical systems are compared to QMG results for physical systems of size 32 x 32 and inverse temperature 
U /T = 128. Comparison of the superfluid density ps and condensed density pc (b) for reference systems of size L = 1 x 1 and 
essentially infinitely large physical systems, cf. Ref. [lol 



where the lattice site indices i have been decomposed 
into an index R, that specifies the cluster and into an in- 
dex a, that specifies the lattice sites within a cluster J^i^ 
Analogously to the physical system, the hoping integrals 
are t'^fj = ~t' and t'^^ = -t' exp[i A (rR^ - rR^)e^] for 
nearest-neighbor hopping processes along the By- and 
the e^^-direction, respectively, and zero otherwise. In 
our calculation, we use the chemical potential /i' and 
the source coupling strength of the reference system 
as variational parameters in the optimization prescrip- 
tion. Since the reference system is complex valued, the 
source coupling strength fa is complex valued too, i.e., 
fa ~ exp[0Q]. Thus, in general, 2L variational pa- 
rameters have to be considered, where L is the number 
of cluster sites. However, for different cluster sites a the 
source coupling strengths fa are interrelated, as can be 
seen from mean field arguments, leading effectively to 
two variational parameters |/| and 0, which we use — in 
addition to the chemical potential /i' — to treat complex 
valued reference systems. 

In Fig. [T] we present the superfluid density ps for differ- 
ent sizes of the reference system ranging from L = 1 x 1 , 
over L = 2xl,toL = 2x2 and essentially infinitely large 
physical systems. For the largest cluster we restrict the 
variational search space to real valued order parameters, 
i.e., we set = 0. Figure [T](a) demonstrates that this 
choice leads to comparable results as obtained with the 
full variational space. Yet, for the restricted variational 
space the computational effort as well as the numerical 
complexity is reduced, since the reference system remains 
real valued. Figure [T](a) shows the superfluid density ps, 
as a function of the chemical potential ^/U evaluated for 
fixed hopping strength t/U — 0.02. The chemical poten- 
tial ranges from p/U = to ^/U = 3. As the hopping 
strength is small, three regions with ps = are present. 



corresponding to the Mott insulating phase. In between 
these regions, we observe a finite superfluid density ps 
indicating the occurrence of the superfluid phase. In ad- 
dition to the VGA results, we show QMC results with er- 
rorbars (barely visible) for physical systems of size 32 x 32 
and inverse temperature U/T = 128. The QMC calcu- 
lations were performed with the ALPS librarj*^ and the 
ALPS applications)^ Particularly, we use the stochas- 
tic series expansion representation of the partition func- 
tion with directed loop updates^^"— where the super- 
fluid density is evaluated via the winding numberJ^iii 
The superfluid density ps obtained from VGA agrees re- 
markable well with the QMC results. Furthermore, VGA 
results are almost independent of the size L of the refer- 
ence system, signaling convergence to the correct results 
even for i = 1 x 1 site clusters. The superfluid den- 
sity Ps is compared to the condensate density pc = (ai) 
in Fig.[T](b), cf. Ref. [13. It can be observed that the 
superfluid density is always larger than the density of 
the Bose-Einstein condensate. However, the difference 
between the two densities is rather small, since a very 
dilute Bose gas is investigated. 

In Fig. [5] we evaluate (a) the superfluid density ps and 
(b) the superfluid fraction ps/n (n is the particle den- 
sity) for fixed chemical potential p/U = 0.4 as a func- 
tion of the hopping strength t/U. The hopping strength 
ranges from t/U = to t/U — 1, which is already 
very deep in the supcrfiuid phase. For p/U = 0.4 the 
phase boundary between the Mott and the superfluid 
phase is located at t/U w 0.06. In the superfluid phase 
close to the phase boundary the superfluid density rises 
quickly from zero developing an almost linear behavior 
for t/J7 > 0.2. In the latter parameter regime the super- 
fluid fraction is larger than 90% signaling that already 
a very large amount of the lattice bosons is superfluid. 
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FIG. 2: (Color online) Superfluid density ps (a) and superfluid 
fraction ps/n (b) ranging deep in the superfluid phase evalu- 
ated for constant chemical potential n/U = 0.4 as a function 
of the hopping strength t/U. Results obtained by means of 
VGA for reference systems of size L = 1 x 1 and essentially 
infinitely large physical systems are compared to QMG results 
for physical systems of size 32 x 32 and inverse temperature 
U/T = 128. 



As emphasized in Ref. [33, a relatively sharp crossover 
from a strongly-correlated superfluid, characterized by a 
superfluid fraction which is well below 1, to a weakly- 
correlated superfluid, where the superfluid fraction is al- 
most 1, can be observed, see Fig.[2](b). In addition to 
the VGA results evaluated for reference systems of size 
L = 1x1 and essentially infinitely large physical systems, 
we show QMC results for physical systems of size 32 x 32 
and inverse temperature U/T = 128, which again exhibit 
perfect agreement. 

In Fig.[3]wc focus on the quantum critical region close 
to the tip of the first Mott lobe, which is the most chal- 
lenging one. In particular, we evaluate the particle den- 
sity n, the condenstate density pc, and superfluid den- 
sity Ps- In the first row we show results for fixed chem- 
ical potential p/U = 0.4 as a function of the hopping 
strength t/U, whereas in the second row we keep hop- 
ping strength fixed at t/U = 0.05 and vary the chemical 



potential p/U. We compare VGA results with QMG and 
mean-field (MF). The most important observation is that 
MF is far off QMC and VGA. For p/U = 0.4 MF predicts 
the phase transition to be at a much smaller value of t/U 
than QMG and VGA. This leads to significant deviations 
in both the density and condensate density as compared 
to QMG and VGA. For fixed t/U = 0.05 MF does not en- 
ter the Mott region and thus does not predict a plateau in 
the density. For both investigated situations (fixed p/U 
and fixed t/U) the results obtained by means of VGA 
and QMG agree quite well. For the QMG simulations we 
used lattices of size 32 x 32 and inverse temperatures of 
U/T = 128. The VGA results are obtained at zero tem- 
perature for clusters of size 1x1 and 2x1, respectively, 
and essentially infinitely large physical systems. In this 
challenging regime small differences between VGA and 
QMG are observable for the condensate density. For the 
reference system sizes considered here, results are almost 
identical. Larger reference systems might still reduce the 
difference between VGA and QMG. However, close to the 
phase transition finite size and finite temperature effects 
might still be important for the QMG results, and thus a 
proper finite size scaling of these data might also reduce 
the discrepancy between the two approaches. Note that 
for fixed hopping t/U = 0.05 there is a very small region 
at p/U 0.85, where it is difficult to numerically deter- 
mine the stationary point of the grand potential. Such 
a region is also present between the first and the second 
and between the second and the third Mott lobe in Fig.[T] 
However, there it is barely visible since the spacing be- 
tween two consecutive p datapoints is larger than this 
gap. This failure appears to be related to the fact that 
two solutions adiabatically connected to two sectors with 
different particle numbers, i. e. the two neighboring Mott 
regions, meet and try to avoid each other. However, we 
want to emphasize that this problem affects only a tiny 
region of the phase diagram. When keeping the chemi- 
cal potential fixed at /i/C/ = 0.4 solutions can be easily 
found for all values of the hopping strength. 

Finally, we want to emphasize that the VGA results are 
obtained with very modest computational effort and that 
excellent agreement with QMG can be observed, even for 
very small reference systems. 



IV. CONCLUSIONS 

In the present work, we extend the self-energy func- 
tional approach to the U{1) symmetry broken, superfluid 
phase of correlated lattice bosons. A crucial point of this 
extension is the identification of a quantity, termed D, 
which is the companion of the self-energy S in the super- 
fluid phase. We also identify the appropriate (nonuni- 
versal) functional fls which is stationary at the physical 
values of the self-energy S and of D. In analogy to the 
self-energy, which is the difference of the interacting and 
non-interacting Green's function, the quantity D is re- 
lated to the difference of the order parameter of the in- 
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FIG. 3; (Color online) Particle density n (left), condensate density pc (middle), and superfluid density ps (right) evaluated 
around the quantum critical region close to the tip of the first Mott lobe. Comparison of the data obtained by means of VGA 
(for essentially infinitely large physical systems and reference systems as stated in the legends), QMC (for physical systems of 
size 32 X 32 and inverse temperatures U /T = 128), and mean-field. The first row (a.*) shows results for fixed chemical potential 
p,/U = 0.4 as a function of the hopping strength t/U , whereas the second row (b.*) shows results for fixed hopping strength 
t/U = 0.05 as a function of the chemical potential p,/U . 



teracting and non-interacting systems. Thus, D is zero 
in the normal phase and for J7 = 0. From these relations 
also follows that both E as well as D vanish in the non- 
interacting case. Importantly, when the functional lis is 
evaluated at the exact values of S and D it corresponds 
to the grand potential of the physical Hamiltonian. To 
evaluate the functional, we proceed as in the original self- 
energy functional approach^i^ and introduce a reference 
system, which is a cluster decomposition of the physical 
system. Importantly, the reference system shares its two- 
particle interaction with the physical system, and can be 
exactly solved by numerical methods. By comparison of 
the functionals, the universal part of l^s, denoted as J^, 
can be eliminated, which allows to evaluate ^Ig exactly on 
the subspace of S and D, spanned by the possible sets 
of reference systems. The results presented are shown 
to be equivalent to the ones obtained by a more heuris- 
tic method, the pseudoparticle approach introduced in 
and thus provide rigorous variational grounds 
for that approach. In addition, the extended self-energy 
functional approach introduced here allows to envision 
more general reference systems, in which bath sites are 
incorporated to enlarge the space of possible self-energies 
E, and possibly bridge over to (Cluster) Dynamical Mean 



Field Theory (DMFT) J^i^^ For future research it would 
be interesting to verify whether in the limit of an infinite 
number of bath sites and for a single correlated site as a 
reference system, our superfluid SFA becomes equivalent 
to DMFT for superfluid bosons^^^ as it is the case in 
the normal phase For a finite number of bath sites this 
is certainly not the case, since the order parameter in the 
reference system differs from the physical one. 

We also presented how the superfluid density can be 
evaluated by means of this extended variational cluster 
approach. To this end we applied a phase twisting field to 
the system. We evaluated the superfluid density for the 
two-dimensional Bose-Hubbard model and compared the 
extended variational cluster approach results with unbi- 
ased quantum Monte Carlo results, yielding remarkable 
agreement. We want to emphasize that the extended 
self-energy functional approach is not only applicable to 
the Bose-Hubbard model but to a large class of lattice 
models, which exhibit a condensed phase. This includes 
experimentally interesting systems such as disordered 
bosons, multicomponent systems (Bose-Bose mixtures 
or Bose-Fermi mixtures) and light matter systems^ii^ 
Strictly speaking, the method cannot treat long-range 
interactions, such as dipolar ones, exactly4^i^ However, 
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the long-range part can be incorporated on a mean-field 
level4^ In principle, the present approach can be applied 
to systems with broken translational invariance as well, 
and, for example, can consider the effect of a confining 
magnetic trap. However, in this case one has to aban- 
don the Fourier transform in the cluster vectors and work 
in real space and, thus, deal with larger matrices and a 
larger number of variational parameters. A convenient, 
numerically less expensive alternative, is to adopt the so- 
called local density approximation!^ 



b. Trace in t and in Matsubara space 
In T space we have 

TtM = /3-Hr / dr M{t,t + 0+) . 
Jo 

The transformation of M{t,t') to Matsubara space is 
defined as 



w„,a;„ e 
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Appendix A: Notation and conventions 



The inverse transformation reads 

MiiJn,Uj'„) = [ dTdT'M{T,T')e"^" 



Combining the equations above, the trace becomes 
1-0 

= r'^trAfK,w„)e'"""^ . 



Matrix notation 



a. General 



In order to simplify our notation we omit time argu- 
ments, whenever this does not cause ambiguities. There- 
fore, two-point functions such as Green's functions, self- 
energies, etc. arc interpreted as matrices in Nambu, or- 
bital, and T space. One-point objects such as A {A) 
are interpreted as column (row) vectors in the same 
space. Matrix-matrix and vector-matrix products are 
understood throughout, whereby internal r variables are 
considered to be integrated over. In addition, the trans- 
posing operator also acts on time variables. Traces 
Tr contain an integral over t and a trace tr over orbital 

indices, i.e., TrM = tr dTM{T,T + 0+), where the 

0+ leads to the well known convergence factor e"^"° in 
Matsubara space. 

(Functional) derivatives with respect to matrices are 
defined "transposed": 



5X 



5M,,{t',t) 



Finally, there are two types of products between row 
(in the form v) and column (m) vectors, depending on 
the order: On the one hand the product vu produces a 
scalar (all indices are summed/integrated over). On the 
other hand, inverting the order, as in uv gives a matrix, 
as indices are "external" and, thus, not summed over. 



c. Logarithm 

There are some subtle points concerning logarithms of 
two-point functions. Although these issues are imma- 
terial for the final result, we prefer to specify them in 
detail. 

The logarithm of G considered as a matrix in the con- 
tinuum variable t is defined up to an infinite constant 
which depends on the the discretization step 5 (see be- 
low). In addition, the trace of the logarithm carried out 
in Matsubara space diverges as well (despite the conver- 
gence factor e'"""^). The usual result presented in the 
literature (see, for instance Ref. |47| ) implicitly assumes 
that an infinite constant has been subtracted. In order 
to avoid these undetermined infinite terms, we subtract 
them explicitly at the outset with the help of the "infinite 
energy" Green's function 

Go,(t, t') = E Goo(w„)e-*'-"(^-^') 

n 

1 



Goo(a;„) = 1 



iUJn — E ' 



where it is understood that the E — > +oo limit is taken 
at the end of the calculation. This choice guarantees, 
for example, that TrlnG/Goo, where G is the Green's 
function in normal (i.e. not Nambu) notation, vanishes 
in the limit /i — ?> — oo, where fi is the chemical potential. 

The Fourier transform defined in App. lA 1 bl allows to 
define the logarithm of G in r space, apart from an in- 
finite multiplicative constant, which originates from the 
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fact that the Fourier transformation is not and cannot be 
normahzed in the continuum limit. In particular, 

[ln(-G)](r,r') = 5^[ln(-G)]K,a.;)e— 

71, n' 



2. Symmetry of Green's functions and other 
two-point functions 

The action in Eq. ([5]) is invariant under the transforma- 
tion Go — ?► (TGqT), where the transposing operator 
also acts on time variables and T is defined in Eq. ([8]) . 
This is due to the fact that 

Air')G^\r',r)A{r) = A{r'frG^\r',r)rAirf 
= A{r)iTG^\r',rfT)A{r'). 

Therefore, we choose Gq to obey the symmetry 

Go = (TG^r) . (Al) 

The same symmetry is obeyed by other two-point func- 
tions, such as the interacting Green's function G, the 
self-energy E, and their inverse. 

In principle, this redundancy renders relations such as 
Eq. (jlSp non invcrtiblc. In order to avoid this, we adopt 



the convention that functional inversions are carried out 
within the subspacc of two-point functions obeying the 
relation Eq. ()A1|) . In addition, we adopt the following 
convention for functional derivatives of an arbitrary func- 
tional S with respect to a two-point function X: 



3. Continuum limit of the functional integral 

In principle, the expression Eq. (jlOp should be under- 
stood such that adjoint fields a are evaluated at a later 
imaginary time t + 6, whereby S is the width of the dis- 
cretization mesh of the interval (0,/3). The continuum 
limit 5-^0 should be taken after having carried out 
the functional integration, see, e.g. Ref. |4^. Taking this 
limit at the outset amounts to neglecting the so-called 
"contribution from infinity" This can be achieved 
by effectively replacing the normal-ordered Hamiltonian 
with a "symmetrically ordered" one, which is suitably 
symmetrized among possible permutation of creation and 
annihilation operatorsi^ In particular, for the noninter- 
acting part, this amounts to replacing the operator ex- 
pression a^a by ^(a^a -I- aa^) = a^a + ^. Therefore, we 
should keep in mind that the grand-potential fl^ corre- 
sponds to such a symmetrized Hamiltonian. 
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